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Abstract 

There are three key factors of a system of coupled oscillators that characterize the interaction 
among them: coupling (how to affect), delay (when to affect) and topology (whom to affect). For 
each of them, the existing work has mainly focused on special cases. With new angles and tools, 
this paper makes progress in relaxing some assumptions of these factors. There are three main 
results in this paper. First, by using results from algebraic graph theory, a sufficient condition 
is obtained which can be used to check equilibrium stability. This condition works for arbitrary 
topology. It generalizes existing results and also leads to a sufficient condition on the coupling 
function with which the system is guaranteed to reach synchronization. Second, it is known that 
identical oscillators with sin() coupling functions are guaranteed to synchronize in phase on a 
complete graph. Using our results, we demonstrate that for many cases certain structures instead 
of exact shape of the coupling function such as symmetry and concavity are the keys for global 
synchronization. Finally, the effect of heterogenous delays is investigated. We develop a new 
framework by constructing a non-delayed phase model that approximates the original one in the 
continuum limit. We further derive how its stability properties depend on the delay distribution. 
In particular, we show that heterogeneity, i.e. wider delay distribution, can help reach in-phase 
synchronization. 



I. INTRODUCTION 



The system of coupled oscillators has been widely studied in different disciplines ranging 
from biology [IHH] and chemistry [HJ [7j to engineering [8j |9] and physics [TQl tH] • The possible 
behavior of such a system can be complex. For example, the intrinsic symmetry of the 
network can produce multiple limit cycles or equilibria with relatively fixed phases (phase- 
locked trajectories) |12j . which in many cases can be stable [13]. Also, the heterogeneity in 
the natural oscillation frequency can lead to incoherence [14J or even chaos [15]. 

One particular interesting question is whether the coupled oscillators will synchronize in 
phase in the long run [T6T420] . Besides its clear theoretical value, it also has rich applications 
in practice. 

In essence, there are three key factors of a system of coupled oscillators that characterize 
the interaction among oscillators: coupling, delay and topology. For each of them, the 
existing work has mainly focused on special crises dbS explained below. In this paper, further 
research will be discussed on each of these three factors: 



Topology (whom to affect, Section III B ): Current results either restrict to complete 

study local stability of topol- 
or introduce dynamic con- 

2^ 



graph or a ring topology for analytical tractability [19] 
ogy independent solutions over time varying graph [2TH23] 

trailers to achieve synchronization for time-varying uniformly connected graphs 
[25] . We develop a graph based sufficient condition which can be used to check equi- 
librium stability for any fixed topology. It also leads to a family of coupling function 
with which the system is guaranteed to reach global phase consensus for arbitrary 
undirected connected graph using only physically meaningful state variables. 



Coupling (how to affect, Section III C ) : The classical Kuramoto model [2] assumes a 
sin() coupling function. Our study hints that certain symmetry and convexity struc- 
tures are enough to guarantee global synchronization. 



Delay (when to affect, Section IV): Existing work generally assumes zero delay among 
oscillators or require them to be bounded up to a constant fraction of the period [26] . 
This is clearly not satisfactory especially if the oscillating frequencies are high. We 
develop a new framework to study unbounded delays by constructing a non-delayed 
phase model that is equivalent to the original one. Using this result, we show that 
wider delay distribution can help reach synchronization. 



In this paper we study weakly coupled oscillators, which can be either pulse-coupled or 
phase-coupled. Although most of the results are presented for phase-coupled oscillators, they 
can be readily extended for pulse-coupled oscillators (see, e.g., [271 EH])- It is worth noticing 
that results in Section |III| are independent of the strength of the coupling and therefore the 
weak coupling assumption is not necessary there. Preliminary versions of this work has been 
presented in [29j and [3D] . 

The paper is organized as follows. We describe pulse-coupled and phase-coupled oscil- 
lator models, as well as their common weak coupling approximation, in Section [IT] Using 
some facts from algebraic graph theory and potential dynamics in Section HI A we present 
the negative cut instability theorem in Section III B 1 to check whether an equilibrium is 
unstable. This then leads to Proposition [T] in Section HI B 2 which identifies a class of 
coupling functions with which the system always synchronizes in phase. It is well known 



2 



that the Kuramoto model produces global synchronization over a complete graph. In Sec- 
tion III C we demonstrate that a large class of coupling functions, in which the Kuramoto 
model is a special case, guarantee the instability of most of the limit cycles in a complete 
graph network. Section IV is devoted to the discussion of the effect of delay. An equivalent 
non-delayed phase model is constructed whose coupling function is the convolution of the 
original coupling function and the delay distribution. Using this approach, it is shown that 
sometimes more heterogeneous delays among oscillators can help reach synchronization. We 
conclude the paper in Section |V| 



II. COUPLED OSCILLATORS 

We consider two different models of coupled oscillators studied in the literature. The 
difference between the models arises in the way that the oscillators interact between each 
other, and their dynamics can be quite different. However, when the interactions are weak 
(weak coupling), both systems behave similarly and share the same approximation. This 
allows us to study them under a common framework. 

Each oscillator is represented by a phase 9i in the unit circle S 1 which in the absence of 
coupling moves with constant speed 9\ = co. Here, S 1 represents the unit circle, or equivalently 
the interval [0, 2tt] with and 2-7T identified (0 = 2n), and uj = ^ denotes the natural 
frequency of the oscillation. 



A. Pulse-coupled Oscillators 

In this model the interaction between oscillators is performed by pulses. An oscillator j 
sends out a pulse whenever it crosses zero (9j = 0). When oscillator i receives a pulse, it will 
change its position from 9{ to 9i + eKij(Qi). The function K{j represents how other oscillators' 
actions affect oscillator i and the scalar e > is a measure of the coupling strength. These 
jumps can be modeled by a Dirac's delta function 5 satisfying 5(t) = Wt ^ 0, 5(0) = +oo, 
and J 5(s)ds = 1. The coupled dynamics is represented by 

9\(t) =ou + ecoJ2 ~ Vij)), (1) 

where rjij > is the propagation delay between oscillators i and j (77^ = 77^), and Mi is 
the set of i's neighbors. The factor of uj in the sum is needed to keep the size of the jump 
within eKij(9i). This is because 9j(t) behaves like ut when crosses zero and therefore the 
jump produced by 5(9j(t)) is of size f 5{9j{t))dt = co> _1 [28] . 

The coupling function be classified based on the qualitative effect it produces in 

the absence of delay. After one period, if the net effect of the mutual jumps brings a pair 
of oscillators closer, we call it attractive coupling. If the oscillators are brought further 
apart, it is considered to be repulsive coupling. The former can be achieved for instance if 
< for ^ G [0, 7r) and Kij(9) > for 9 £ [tc, 2tt). See Figure [l] for an illustration of an 
attractive coupling Kij and its effect on the relative phases. 

This pulse-like interaction between oscillators was first introduced by Peskin [2] in 1975 
as a model of the pacemaker cells of the heart, although the canonic form did not appear in 
the literature until 1999 [28]. In general, when the number of oscillators is large, there are 
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FIG. 1. Pulse-coupled oscillators with attractive coupling. 



several different limit cycles besides the in-phase synchronization and many of them can be 
stable [T3] . 

The question of whether this system can collectively achieve in-phase synchronization was 
answered for the complete graph case and zero delay by Mirollo and Strogatz in 1990 [20J. 
They showed if «y (9) is strictly increasing on (0, 2ir) with a discontinuity in (which resem- 
bles attractive coupling), then for almost every initial condition, the system can synchronize 
in phase in the long run. 

The two main assumptions of [20J are all to all comunication and zero delay. Whether 
in-phase synchronization can be achieved for arbitrary graphs has been an open problem for 
more than twenty years. On the other hand, when delay among oscillators is introduced the 
analysis becomes intractable. Even for the case of two oscillators, the number of possibilities 
to be considered is large [3H [32] . 

B. Phase-coupled Oscillators 

In the model of phase-coupled oscillators, the interaction between neighboring oscillators 
i and j G H% is modeled by change of the oscillating speeds. Although in general the speed 
change can be a function of both phases (9i, 9j), we concentrate on the case where the speed 
change is a function of the phase differences fij{4>j{t — T]ij) — fait)). Thus, since the net speed 
change of oscillator i amounts to the sum of the effects of its neighbors, the full dynamics is 
described by 

Ut) = u + e ~ Vij) - &(*))■ (2) 

The function is usually called coupling function, and as before rjij represents delay and 
Mi is the set of neighbors of i. 
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FIG. 2. Phase-coupled oscillators with attractive and repulsive coupling. 



A similar definition for attractive and repulsive couplings can be done in this model. We 
say that the coupling function is attractive if, without delays, the change in speeds brings 
oscillators closer, and repulsive if they are brought apart. Figure [2] shows typical attractive 
and repulsive coupling functions where arrows represent the speed change produced by the 
other oscillator; if the pointing direction is counter clockwise, the oscillator speeds up, and 
otherwise it slows down. 

When fij = |vsin(), K > (attractive coupling), this model is known as the classical 
Kuramoto model [33]. Intensive research work has been conducted on this model, however 
convergence results are usually limited to cases with all to all coupling (Mi = Af\{i}, i.e., 
complete graph topology) and no delay (rjij = 0), see e.g. [T9l [33] , or to some regions of the 
state space [26] . 



C. Weak Coupling Approximation 

We now concentrate in the regime in which the coupling strength of both models is weak, 
i.e. 1 ^> e > 0. For pulse-coupled oscillators, this implies that the effect of the jumps 
originated by each neighbor can be approximated by their average [27J. For phase-coupled 
oscillators, it implies that to the first order <pi(t — is well approximated by <p%(t) — curjij. 

The effect of these approximations allows us to completely capture the behavior of both 
systems using the following equation where we assume that every oscillator has the same 
natural frequency u and only keep track of the relative difference using 

<j>i = fij{<f>j - <t>i- (3) 
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For pulse-coupled oscillators, the coupling function is given by 



fM = ^H), (4) 

and the phase lag ipij = urjij represents the distance that oscillator i's phase can travel 
along the unit circle during the delay time r]ij. Equation Q also shows that the attrac- 
tive/repulsive coupling classification of both models are in fact equivalent, since in order to 
produce the same effect Kij and should be mirrored, as illustrated in Figure [T] and Figure 

EJ 

Equation ^ captures the relative change of the phases and therefore any solution to ^ 
can be immediately translated to either Q or ^ by adding cot. For example, if (ft* is an 
equilibrium of ([3]), by adding cot, we obtain a limit cycle in the previous models. Besides 
the delay interpretation for ipij, (|3j) is also known as a system of coupled oscillators with 
frustration, see e.g. [35] . 

From now on we will concentrate on (|3j) with the understanding that any convergence 
result derived will be immediately true for the original models in the weak coupling limit. 
We are interested in the attracting properties of phase-locked invariant orbits within T , 
which can be represented by (p(t) = co*Hn + (ft*, where ljy = (1, • • • , 1) T G 7~ N , and 0* and 
co* are solutions to 

&($-#-^-),vz. ( 5 ) 

Whenever the system reaches one of these orbits, we say that it is synchronized or phase- 
locked. If furthermore, all the elements of <fi* are equal, we say the system is synchronized 
in-phase or that it is in-phase locked. It is easy to check that for a given equilibrium cff of 
(J3|, any solution of the form <p* + AIat, with A G R, is also an equilibrium that identifies the 
same limit cycle. Therefore, two equilibria <fi l >* and <fi 2 '* will be considered to be equivalent, 
if both identifies the same orbit, or equivalently, if both belongs to the same connected set 
of equilibria 

E r :={</» eT N \<j> = 4>* + AIat, Ag R}. (6) 



III. EFFECT OF TOPOLOGY AND COUPLING 

In this section we concentrate on the class of coupling function that are symmetric 
(fij — fji Vij), odd ( fij(—0) = —fij(0)) and continuously differentiable. We also assume 
that there is no delay within the network, i.e. ipy = Vij. Thus, ^ reduces to 

k = fiMj - <Pi)- (7) 

In the rest of this section we progressively show how with some extra conditions on we 
can guarantee in-phase synchronization for arbitrary undirected graphs. Since we know that 
the network can have many other phase-locked trajectories besides the in-phase one, our 
target is an almost global stability result [36], meaning that the set of initial conditions 
that does not eventually lock in-phase has zero measure. Latter we show how most of the 
phase-locked solution that appear on a complete graph are unstable under some general 
conditions on the structure of the coupling function. 
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A. Preliminaries 



We now introduce some prerequisites used in our later analysis. 



1. Algebraic Graph Theory 

We start by reviewing basic definitions and properties from graph theory J37J EH] that are 
used in the paper. Let G be the connectivity graph that describes the coupling configuration. 
We V{G) and E(G) to denote the set of vertices (i or j) and undirected edges (e) of G. An 
undirected graph G can be directed by giving a specific orientation a to the elements in the 
set E{G). That is, for any given edge e G E(G), we designate one of the vertices to be the 
head and the other to be the tail giving G° ' . 

Although in the definitions that follow we need to give the graph G a given orientation a, 
the underlying connectivity graph of the system is assumed to be undirected. This is not 
a problem as the properties used in this paper are independent of a particular orientation a 
and therefore are properties of the undirected graph G. Thus, to simplify notation we drop 
the superscript a from G a with the understanding that G is now an induced directed graph 
with some fixed, but arbitrarily chosen, orientation. 

We use P = (V~,V + ) to denote a partition of the vertex set V(G) such that V(G) = 
V~ U V + and V~ fl V + = 0. The cut C{P) associated with P, or equivalently C(V~, V + ), 
is defined as C(P) := {ij G E(G)\i G V~,j G V + , or vice versa.}. Each partition can be 
associated with a vector column cp where cp(e) = 1 if e goes form V" to V + , cp(e) = — 1 
if e goes form V + to V~ and cp(e) = if e stays within either set. 

There are several matrices associated with the oriented graph G that embed information 
about its topology. However, the one with most significance to this work is the oriented 
incidence matrix B G Rl^( G )l x l E ( G )l where B(i,e) — 1 if i is the head of e, B(i,e) = —1 if i 
is the tail of e and B(i,e) = otherwise. 



2. Potential Dynamics 

We now describe how our assumptions on /j,- not only simplifies considerably the dynam- 



ics, but also allows us to use the graph theory properties introduced in Section |III A 1| to 
gain a deeper understanding of 

While fij being continuously different iable is standard in order to study local stability and 
sufficient to apply LaSalle's invariance principle [39j, the symmetry and odd assumptions 
have a stronger effect on the dynamics. 

For example, under these assumptions the system ([7]) can be compactly rewritten in a 
vector form as 

= —eBF(B T 4>) (8) 



where B is the adjacency matrix defined in Section III A 1 and the map F : S(G) — > S(G) is 

E{y) = {fij{Vij))ij<EE{G)- 

This new representation has several properties. First, from the properties of B one can 
easily show that ^ can only hold with to* = for arbitrary graphs [16] (since Nui* = 
o;*1^-1jv = —sl^jBF{B T (f)) = 0), which implies that every phase-locked solution is an 
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equilibrium of ([7b and that every limit cycle of the original system ^ can be represented 
by some E\ on |7|. 

However, the most interesting consequence of (j8() comes from interpreting F(y) as the 
gradient of a potential function 



w (y)= E r Us)ds. 

JO 



ijeE(G) 

Then, by defining V(0) = (W o £? T )(0) = W(B T (j)), (J8| becomes a gradient descent law for 
V{<f>), i.e., 

= -sBF(B T (j)) = -sBVW(B T <f)) = -eW(0), 

where in the last step above we used the property V(W o B T )((f)) = BVW \B T 0) . This 
makes V(<f>) a natural Lyapunov function candidate since 

V(<f>) = (W(0),0> = -s |W(0)| 2 = — |0| 2 < 0. (9) 



Furthermore, since the trajectories of (|8j) are constrained into the iV-dimensional torus 
T , which is compact, V(<f>) satisfies the hipotesis of LaSalle's invariance principle (Theorem 
4.4 [39J), i.e. there is a compact positively invariant set. T A and a function V : T A — > R that 
decreases along the trajectories 4>(t). Therefore, for every initial condition, the trajectory 
converges to the largest invariant set M within {V = 0} which is the equilibria set E = 

{^er N \<p^0} = [j^E r . 



Remark 1. The fact that symmetric and odd coupling induces potential dynamics is well 
know in the physics community J^0| /. However, it has been also rediscovered in the control 
community [H\[ for the specific case of sine coupling. Clearly, this is not enough to show 
almost global stability, since it is possible to have other stable phase-locked equilibrium sets 
besides the in-phase one. However, if we are able show that all the non-in-phase equilibria 
are unstable, then almost global stability follows. That is the focus of the next section. 



B. Negative Cut Instability Condition 



We now present the main results of this section. Our technique can be viewed as a 
generalization of [19]. By means of algebraic graph theory, we provide a better stability 
analysis of the equilibria under a more general framework. We further use the new stability 
results to characterize that guarantees almost global stability. 



1. Local Stability Analysis 



In this section we develop the graph theory based tools to characterize the stability of each 
equilibrium. We will show that given an equilibrium cp* of the system ([8]), with connectivity 
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graph G and /y as described in this section. If there exists a cut C(P) such that the sum 



£ 4 



>*) < o, 



(10) 



iieC(P) 



the equilibrium 0* is unstable. 

Consider first an equilibrium point 0*. Then, the first order approximation of ^ around 
6* is 



-sB 



d_ 
dy 



F(B J 



B 1 



were 



— 0* is the incremental phase variable, and ^F(B T c l 



Jacobian of F(y) evaluated at B T 
Now let A = -sB 



i.e. 



9 -F(B^ 



*) g R l^(G)|x|£(G)| is the 

'j ~ #)}tfeJS(G)) • 



it is possible to numerically calcu 



By- v~ r , = diag ({/£ 
_B T and consider the linear system 50 = A5<f). Although 

ate the eigenvalues of A given 0* to study the stability, 
here we use the special structure of A to provide a sufficient condition for instability that 
has nice graph theoretical interpretations. 

Since A is symmetric, it is straight forward to check that A has at least one positive 
eigenvalue, i.e. 0* is unstable, if and only if x 1 Ax > 0. Now, given any partition P = 
(V~,V + ), consider the associated vector cp, define xp such that Xj = \ if % G V + and 
Xi = — | if i G V' 
that 



Then it follow from the definition of B that cp = B xp which implies 



-1 



x T P Axp 



c P 



d_ 
dy 



F(B T <f) 



c P 



ijeC(P) 



Therefore, when condition (10) holds, A 
real part is positive. 



-eBDB T has at least one eigenvalue whose 



Remark 2. Equation (10) provides a sufficient condition for instability; it is not clear 



what happens when (10) does not hold. However, it gives a graph- theoretical interpretation 



that can be used to provide stability results for general topologies. That is, if the minimum 
cut cost is negative, the equilibrium is unstable. 



Remark 3. Since the weights of the graph fij{4>*j —<P*) are functions of the phase difference, 
(10) holds for any equilibria of the form 0* + Aljv- Thus, the result holds for the whole set 
defined in 

When (10) is specialized to P = ({i}, V(G)\{i}) and fij{0) = sin(0), it reduces to the 
instability condition in Lemma 2.3 of |19j : i.e., 



< 0. 



(11) 



However, (10) has a broader applicability spectrum as the following example shows. 



Example 1. Consider a six oscillators network as in Figure^ where each node is linked 
with its four closest neighbors and fij{9) = sm(9). Then, by symmetry, it is easy to verify 
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that 




FIG. 3. The network of six oscillators (Example 4) 



7T 




FIG. 4. Unstable equilibrium eft*. Initial condition cpQ = (j)* + 



We first study the stability of (J)* using (11) as in JTPj /. By substituting (12) in cos(0* — 0*) 
Vzj G -B(G) we /md t/iat t/ie edge weights can only take two values: 



COS ( 

cos( 



3/ 
3 > 



1 

2> 



1 

' 2 • 



i/ j = i ± 1 mod 6 
if j = i ± 2 mod 6 



T/ien, since am/ cut t/iat isolates one node from the rest (like C\ = C({1}, V(G)\{1}) in 
Figure^ will always have two edges of each type, their sum is zero. Therefore, ( [TTj ) cannot 
be used to determine stability. 
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If we now use condition (10) instead, we are allowed to explore a wider variety of 
cuts that can potentially have smaller costs. In fact, if instead of C\ we sum over C2 = 
C({1,2,6},{3,4,5}) ; we obtain, 



COS 



-1 < 0, 



which implies that (f>* is unstable. 

Figure ^\ verifies the equilibrium instability. By starting with an initial condition O — 
0* + 5<p close to the equilibrium <f>* , we can see how the system slowly starts to move away 
from <f>* towards a stable equilibrium set. 

Furthermore, we can study the whole family of non-isolated equilibria given by 



n 2n Aix 5tv 

£1, g + £2, -y + £3, TT + £1, — + £2, y + £3 



1 T 



(13) 



where £i,£2, £3 G R, which due to Remark^ we can reduce (13) to 



7r 2ty 4f 5f 

0, - + M, — + A 2 ,vr, — + Ai, — + A 2 



(14) 



with Ai = £2 — £1 and A2 = £3 — £1. 




FIG. 5. Minimum cut value C*(Ai,A2) showing that the equilibria (13) are unstable 



Instead of focusing on only one cut, here we compute the minimum cut value (10) over 
the 31 possible cuts, i.e. C*(Ai,A 2 ) := mirip Y^ijec(p) flji^ji^ij ^2)* — 0,*(Ai,A 2 )). Figure 
[| show the value of C*(\i, \ 2 ) for Aj G [— 7r,7r]. Since C*(Ai,A 2 ) is 27T-periodic on each 
variable and its value is negative for every Ai,A 2 G [— 7T,7r], the family of equilibria (14) 



(and consequently {13)) is unstable 
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2. Almost Global Stability 



Condition (10) also provides insight on which class of coupling functions can potentially 
give us almost global convergence to the in-phase equilibrium set E% N - If it is possible to 
find some with /^(0) > 0, such that for any non-in-phase equilibrium 0*, there is a cut 
C with Yliij&c fiji&j ~~ $) < 0' then the in-phase equilibrium set will be almost globally 
stable |13j . The main difficulty is that for general and arbitrary network G, it is not easy 
to locate every phase-locked equilibria and thus, it is not simple to know in what region of 
the domain of /y the slope should be negative. 

We now concentrate on the one-parameter family of functions J-^, with b G (0, 7r), such 
that fij G Tb whenever /y is symmetric, odd, continuously differentiable and 

• /£.(0; b) > 0, W G (0, b) U (2tt - b, 2vr), and 

• fffib) <0, W9e (6, 27T-6). 

See Figure [2] for an illustration with b = | . Also note that this definition implies that if 
6) G the coupling is attractive and fij (9; b) > \/9 G (0, tt). This last property will 
be used later. We also assume the graph G to be connected. 

In order to obtain almost global stability we need b to be small. However, since the 
equilibria position is not known a priori, it is not clear how small b should be or if there is 
any b > such that all nontrivial equilibria are unstable. We therefore first need to estimate 
the region of the state space that contains every non-trivial phase-locked solution. 

Let I be a compact connected subset of S 1 and let 1(1) be its length, e.g., if / = S 1 then 
1(1) = 2ir. For any S C V(G) and G T N , define d(4>, S) as the length of the smallest 
interval / such that 0j G I Vz G S, i.e. 

d(<f>,S) = l(I*)= min 1(1). 

Using this metric, together with the aid of Proposition 2.6 of [16] we can identify two 
very insightful properties of the family J-j, whenever the graph G is connected. 



Claim 1. If 0* is an equilibrium point of ([8]) with d(<p* ,V(G)) < n, then either <fi* is an 
in-phase equilibrium, i.e. (f>* = XI n for X G R, or has a cut C with — 0*) < G C. 

Proof. Since d((p*,V(G)) < it, all the phases are contained in a half circle and for the 
oscillator with smallest phase io, all the phase differences (0* — 0* o ) G [0, 71*]. However, 
since f%j(-',b) G T\, implies fij(9;b) > W9 G [0, 7r] with equality only for 9 G {0,7r}, 
^ = E j£ m - = can only hold if 0* - 0*, G {0,tt} Vj G M- Now let 

r = {i£ V^(G?) : rf(0*,{z,z o }) = 0} and V + = V(G)\V~. If V~ = V(G), then 0* is an 
in-phase equilibrium. Otherwise, Vij G C(V _ , V^ + ), fij(4>* - 0*) = /»j(7r) < 0. □ 

We are now ready to establish a bound on the value of b that guarantees the instability 
of the non-in-phase equilibria. 

Claim 2. Consider fij(-; b) G Vij G E(G) and arbitrary connected graph G. Then for any 
b < -^-j- and non-in-phase equilibrium 0*, there is a cut C with fij(4>*j — 0*; b) < 0, Vij G C 
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Proof. Suppose there is a non-in-phase equilibrium 0* for which no such cut C exists. Let 
Vq = {i } and Vq = V(G)\{i } be a partition of V{G) for some arbitrary node i . 

Since such C does not exists, there exists some edge i ji G C(Vq , Vq~), with j\ G V^ + , 
such that — ^io'^) — 0- Move ji from one side to the other of the partition by 

defining Vf := 7 " U {jj and := Now since /^(^ - 0* o ; 6) > 0, then 

In other words, both phases should be within a distance smaller than b. 

Now repeat the argument k times. At the k th iteration, given V k T_ 1 , V fc + 1; again we can find 
some i k -i G V^, J k e V+_ t such that £ fc _u"fc e and fl h _ dk (<J>* jk - <P* k ^b) > 0. 

Also, since at each step d(<p*, {ik-i,jk}) < 

Thus by solving the recursion we get: d(<f>*, V^) < kb. 

After N - 1 iterations we have V^_ x = V(G) and d(<p\V(G)) < (N - 1)6. Therefore, 
since b < jtzt, we obtain 

d(<P*,V(G))<(N-l) J ^ = n. 

Then, by Claim[T]0* is either an in-phase equilibrium or there is a cut C with flj{<j>j — <f>*) < 
Vij G C. Either case gives a contradiction to assuming that 0* is a non-in-phase equilibrium 
and C does not exists. Therefore, for any non-in-phase <fi* and b < jrzr, we can always find 
a cut C with ^(0* - 0*; 6) < 0, G C. □ 



Claim pi allows us to use our cut condition (10) on every non-in-phase equilibrium. Thus, 



since rt8j) is a potential dynamics (c.f. Section III A 2 ), from every initial condition the system 



converges to the set of equilibria E. But when b < -^—^ the only stable equilibrium set inside 
E is the in-phase set E\ N . Thus, E 1n set is globally asymptotically stable. We summarized 
this result in the following Proposition. 

Proposition 1 (Almost global stability). Consider fij(0;b) G and an arbitrary con- 
nected graph G. Then, if b < jr^i, the in-phase equilibrium set E 1n is almost globally 
asymptotically stable. 

This result provides a sufficient condition for almost global asymptotic stability to the in- 
phase equilibrium set E 1n . Although found independently, the same condition was proposed 
for a specific piecewise linear in [H]. Here we extend [JT] in many aspects. For example, 
instead of assuming equal coupling for every edge, our condition describes a large family of 
coupling functions Fb where each f^ can be taken independently from Tb- Also, in [H] the 
construction of fij{0) assumes a discontinuity on the derivative at 6 = b. This can pose a 
problem if the equilibrium 0* happens to have phase differences 0* — 0* = b. Here we do 
not have such problem as is continuously differentiable. 

The condition b < jr-r implies that, when iV is large, should be decreasing in most of 
it domain. Using (pi) this implies that «y should be increasing within the region (6, 2n — b), 
which is similar to the condition on [20] and equivalent when b — > 0. Thus, Proposition [T] 
confirms the conjecture of [20] by extending their result to arbitrary topologies and a more 
realistic continuous for the system ([!]) in the weak coupling limit. 
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C. Complete Graph Topology with a Class of Coupling Functions 



In this subsection we investigate how conservative the value of b found in Section III B 2 
is for the complete graph topology. We are motivated by the results of [19J where it is 
shown that f{6) = sin(6>) (b = |) with complete graph topology ensures almost global 
synchronization. 

Since for general / it is not easy to characterize all the possible equilibria of the system, 
we study the stability of the equilibria that appear due to the equivalence of ^ with respect 
to the action group Sn x T 1 , where Sn is the group of permutations of the N coordinates 
and T 1 = [0,27r) represents the group action of phase shift of all the coordinates, i.e. the 
action of 5 G T 1 is fa (->• fa + 5 Wi. We refer the readers to [12] and [TH] for a detailed study 
of the effect of this property. 

These equilibria are characterized by the isotropy subgroups T of Sn x T 1 that keep them 
fixed, i.e., 70* = fa V7 G T. In [12J it was shown that this isotropy subgroup takes the form 
of 

{S ko x S kl x ••• x S k ) m x 



where fcj and m are positive integers such that (k + k± + ■ ■ ■ + k\ B -\)m = N, Sj is the 
permutation subgroup of Sn of j'-many coordinates and Z m is the cyclic group with action 
fa | — > fa + — ■ The semiproduct x represents the fact that Z m does not commute with the 
other subgroups. 



In other words, each equilibria with isotropy (Sk x S kl x • • • x S k , - ) m x Z m is conformed 
by Ib shifted constellations C\ (I G {0, 1, ... Ib — 1}) of m evenly distributed blocks, with k\ 
oscillators per block. We use 5i to denote the phase shift between constellation Co and C\. 
See Figure [6] for examples these types of equilibria. 




FIG. 6. Equilibria with isotropy (Sk x S kl x Sk 2 ) x -^4 (left) an d (Sk) x Z$ (right) 



Here we will show that under mild assumptions on / and for b — | most of the equilibria 
found with these characteristics are unstable. We first study all the equilibria with m even. 
In this case there is a special property that can be exploited. 
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That is, when / G such that / is even around f , we have 



9 m {8):=JZf(—j + S) (15) 
*— ' m 

3=0 

W2-1 9 

i=o 

m/2-l 

/(-i + fl + /(( T + -J + *)- 5 ) 

i=o 

m/2-l 

,27T. ,2tt 



E /(-j +*)+/(-&•+*)) 

3=0 

E f(-j + s)-f(-j + s) = o 

^— ' m m 



m m 

3=0 

where the third step comes from / being even around tt/2 and 27r-periodic, and the fourth 
from / being odd. 

Having g m {8) = is the key to prove the instability of every equilibria with even m. 
It essentially states that the aggregate effect of one constellation C\ on any oscillator j G 
V(G)\Ci is zero when m is even, and therefore any perturbation that maintains Ci has null 
effect on j. This is shown in the next proposition. 




Proposition 2 (Instability for even m ). Given an equilibrium <fi* with isotropy (S^ X Sk 2 x 
• • • X jSfc, ) m x Z m and / G Til even around |. Then, if m is even, 0* zs unstable. 



Proof. We will show the instability of 0* by finding a cut of the network satisfying (10). 
Let Vq C V(G) be the set of nodes within one of the blocks of the constellation Co and 
consider the partition induced by V , i.e. P = (Vq, V(G)\Vq). Due to the structure of </>*, 
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(10) becomes 



i B 

£ m - 4>*) = -hf'(o) + £ k ^' m (Si), 

ijeC(P) 1=1 

where g' m (S) is the derivative of g m and 5i is the phase shift between the C and C\. Finally, 
since by assumptions g m (d~) = \/S then it follows that g' m (S) = and 

£ &(#-«) = -*i/'(0)<0. 

iieC(p) 



Therefore, by (10), 0* is unstable. □ 

The natural question that arises is whether similar results can be obtained for m odd. 
The main difficulty in this case is that g m {$) = does not hold since we no longer evaluate 
/ at points with phase difference equal to tc such that they cancel each other. Therefore, an 
extra monotonicity condition needs to be added in order to partially answer this question. 
These conditions and their effects are summarized in the following claims. 

Claim 3 (Monotonicity). Given f G Ttl such that f is strictly concave for 9 G [0, 7r], then 

f'(0) - f'(0 -<f>)< 0, <6 -0 <6 <n (16) 
f{9) - f{9 + 0)<O, -7T<6<6 + <f)<0 (17) 

Proof. The proof is a direct consequence of the strict concavity of /. Since f(8) is strictly 
concave then basic convex analysis shows that f'(9) is strictly decreasing within [0, 71*]. There- 
fore, the inequality (16) follows directly from the fact that 9 G [0, tt],9 — (j) G [0, ir] and 
9 — (f) < 9. To show (17) it is enough to notice that since / is odd ( / G J-n), f is strictly 
convex in [it, 2tt]. The rest of the proof is analogous to (16). □ 

Claim 4 (/' Concavity). Given f G J-i such that f is strictly concave for 9 G 
Then for all m > A, f'(^) > |/'(0). 

Proof. Since f'{9) is concave for 9 G [— 7r,7r] then it follows 

f(-) = /'(A m + (1 - \ m f~) > X m f (0) + (1 - \a)f'Q) > Xnf (0) 

where A m = m ^ L - Thus, for m > 4, A m > | and 

f(-) > J/'(0) 

as desired. □ 



7T 7Tl 
2 ' 2 J 



Now we show the instability of any equilibria with isotropy (S^ x x • • • x ) m x Z„ 
for m odd and greater or equal to 7. 



Proposition 3 (Instability for m > 7 and odd). Suppose f G with f concave in [0, ir] 
and f concave in [— f , f ], i/ien /or a// m = 2A; + 1 m& > 3 t/ie equilibria 0*s wzt/i isotropy 
(S^ x ^2 x • • • x "Sjfej ) m x Z m are unstable. 
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FIG. 8. Cut used in Proposition [3} The dots in red represent all the oscillators of some maximal 
set S with d{(j)*,S) < % 

The proof of Proposition [3] also uses our cut condition to show instability, but with a 
different cut induced by the partition P = (S,V(G)\S) of V{G) where S is set to the a 
maximal subset of V{G) such that g?(0, S) < see Figure || for an illustration of P. Notice 
that any of these partitions will include all the oscillators oftwo consecutive blocks of every 
constellation. The details of the proofs are rather technical and are relegated to Appendix 



IV. EFFECT OF DELAY 

Once delay is introduced to the system of coupled oscillators, the problem becomes fun- 
damentally harder. For example, for pulse-coupled oscillators, the reception of a pulse no 
longer gives accurate information about the relative phase difference A<pij = <pj — 0j between 
the two interacting oscillators. Before, at the exact moment when i received a pulse from 
j, <pj was zero and the phase difference was estimated locally by i as A&j = — 0j. How- 
ever, now when i receives the pulse, the difference becomes A0jj = — 0j — ip^. Therefore, 
the delay propagation acts as an error introduced to the phase difference measurement and 
unless some information is known about this error, it is not possible to predict the behavior. 
Moreover, as we will see later, slight changes in the distribution can produce nonintuitive 
behaviors. 

Even though it may not be satisfactory for some applications, many existing works choose 
to ignore delay, (see for e.g., [T8-20J ) . That is mainly for analytical tractability. On the 
other hand, when delay is included [26] the studies concentrate on finding bounds on delay 
that maintain stability. 

In this section we study how delay can change the stability in a network of weakly coupled 
oscillators. A new framework to study these systems with delay will be set up by constructing 
an equivalent non-delayed system that has the same behavior as the original one in the 
continuum limit. We then further use this result to show that large heterogeneous delay 
can help reach synchronization, which is a bit counterintuitive and significantly generalizes 
previous related studies [281 H21 H2] . We will assume complete graph to simplify notation 
and exposition although the results can be extended for a boarder class of densely connected 
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networks. 

The contribution of this section is two fold. First, it improves the understanding of the 
effect of delays in networks of coupled oscillators. And second, it opens new possibilities of 
using delay based mechanisms to increase the region of attraction of the in-phase equilib- 
rium set. We shall build on existing arguments such as mean field approximation |33j and 
Lyapunov stability theory [UJ [19] while looking at the problem from a different perspective. 



A. Mean Field Approximation 

Consider the case where the coupling between oscillators is all to all and identical {Mi = 
Af\{i}, Vi G J\f and fij = f Vz,j). And assume the phase lags ipij are randomly and 
independently chosen from the same distribution with probability density g{ip). By letting 
N +00 and e — > while keeping eN =: e a constant, ^ becomes 

/7T /* + 00 
/ f(a - - ip)g(ip) P (a, t)d^da, (18) 
-7T JO 

where t) is a time- variant normalized phase distribution that keeps track of the fraction 
of oscillators with phase <fi at time t, and v(<p,t) is the velocity field that expresses the net 
force that the whole population applies to a given oscillator with phase <fi at time t. Since 
the number of oscillators is preserved at any time, the evolution of t) is governed by the 
continuity equation 



with the boundary conditions p(0,t) = p(2n,t). Equations ( 18 )-( 19 ) are not analytically 
solvable in general. Here we propose a new perspective that is inspired by the following 
observation. 

Consider the non-delayed system of the form 

fa = u + e'%2H(<f> j -<f> i ), (20) 

where 

r+00 

H{6) = f * g{9) = f{6- i/>)gty)dil> (21) 
Jo 

is the convolution between / and g. 

By the same reasoning of (18) it is easy to see that the limiting velocity field of (20) is 

r-2- 



Vh(4>, t) = uj + e / H(cr — 4>)p{a 1 t)da 
Jo 



'0 

p2n / r+00 \ 

U+ Jo \Jo ~ ^ ~ ^)9^) d ^j Pfo t)dxr 

uj + e / f(a - (j)-ijj)g{ijj)p{a,t)dijjda 



Jo 
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where in the first and the third steps we used (21) and (18) respectively. Therefore, (J3T) and 



(20) have the same continuum limit. 



Remark 4. Although (20) is quite different fro m p| , both systems behave exactly the same 
in the continuum limit. Therefore, as N grows, (20) starts to become a good approximation 
°f PI and therefore can be analyzed to understand the behavior of p|. 




FIG. 9. Effect of delay in coupling shape 



Figure [9] shows how, the underlying delay (in this case the delay distribution) determines 
what type of coupling (attractive or repulsive) produces synchronization. The original func- 
tion / produces repulsive coupling, whereas the corresponding H is attractive. In fact, as we 
will soon see, the distribution of delay not only can qualitatively affect the type of coupling 
but also can change the stability of certain phase-locked limit cycles. 

We now study two example to illustrate how this new approximation can provide signif- 
icant information about performance and stability of the original system. We also provide 
numerical simulations to verify our predictions. 

B. Kuramoto Oscillators 

We start by studying an example in the literature jH] to demonstrate how we can use the 
previous equivalent non-delayed formulation to provide a better understanding of systems 
of coupled oscillators with delay. When f(8) = Ksin(9), H(9) can be easily calculated: 

r+oo 

H(8)= Ksm(6 -ip)g(ip)dip 
Jo 

r+oo r+oo 

= K 3[e i(0 -^£(V)]# = K^[e w / e~ ilp g(^)d^} 
Jo Jo 

= K% [e ie Ce-^] = KCsm(6 - f) 

where 3 is the imaginary part of a complex number, i.e. 3 [a + ib] = b. The values of C > 
and £ are calculated using the identity 

r+oo 

Ce^= / e^g(tp)dip. 
Jo 

This complex number, usually called "order parameter", provides a measure of how the 
phase-lags are distributed within the unit circle. It can also be interpreted as the center of 
mass of the lags ipifs when they are thought of as points (e 4 ^* 3 ') within the unit circle S 1 . 
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Thus, when C ~ 1, the ipi/s are mostly concentrated around £. When C 
distributed such that e 1 ^' rs 0. 



In this example, (20) becomes 



+ e#C ^ sin( 



0- 



0, the delay is 



(22) 



Here we see how the distribution of g{ip) has a direct effect on the dynamics. For ex- 
ample, when the delays are heterogeneous enough such that C ~ 0, the coupling term dis- 
appears and therefore makes synchronization impossible. A complete study of the system 
under the context of superconducting Josephson arrays was performed [H] for the complete 
graph topology. There the authors characterized the condition for in-phase synchroniza- 
tion in terms of K and Ce^. More precisely, when KCe 1 ^ is on the right half of the plane 
(KC cos(^) > 0), the system almost always synchronizes. However, when KCe 1 ^ is on the 
left half of the plane (KC cos(£) < 0), the system moves towards an incoherent state where 
all of the oscillators' phases spread around the unit circle such that its order parameter, i.e. 



N 



becomes zero. 
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FIG. 10. Delay distributions and their order parameter Ce 1 ^ 



We now provide simulation results to illustrate how (22) becomes a good approximation 



of the original system when N is large enough. We simulate the original repulsive (K < 0) 
sine-coupled system with heterogeneous delays and its corresponding approximation (22). 



Two different delay distributions, depicted in Figure [TOj were selected such that their cor- 
responding order parameter lie in different half-planes. 

The same simulation is repeated for N = 5, 10, 50. Figure [TT] shows that when N is 
small, the phases' order parameter of the original system (in red/blue) draw a trajectory 
which is completely different with respect to its approximation (in green). However, as N 
grows, in both cases the trajectories become closer and closer. Since K < 0, the trajectory 
of the system with wider distribution (Ccos£ < 0) drives the order parameter towards the 
boundary of the circle, i.e., heterogeneous delay leads to homogeneous phase. 



C. Effect of Heterogeneity 



We now explain a more subtle effect that heterogeneity can produce. Consider the system 



in (20) where H odd and continuously differentiate. Then, from Section III, all of the 
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FIG. 11. Repulsive sine coupling with heterogeneous delays 



oscillators eventually end up running at the same speed ui with fixed phase difference such 
H{(j)^ — (j)j) cancels Wi. Moreover, we can apply (10) to assess the stability 

H'(<b*-(b*) < 



that the sum ^\ 

of these orbits. Therefore, if we can find a cut C of the network such that J2ijec 
0, the phase-locked solution will be unstable. 

Although this condition is for non-delayed phase-coupled oscillators, the result of this 
section allows us to translate it for systems with delay. Since H is the convolution of the 
coupling function / and the delay distribution function g, we can obtain H'((p* — <$*) < 0, 
even when f'(4>* — <$*) > 0. This usually occurs when the convolution widens the region 
with a negative slope of H. See Figure [9] for an illustration of this phenomenon. 



Figures 12 and 13 show two simulation setups of 45 oscillators pulse-coupled all to all. 
The initial state is close to a phase locked configuration formed of three equidistant clusters 
of 15 oscillators each. The shape of the coupling function / and the phase lags distributions 
are shown in part a. We used Q to implement the corresponding pulse-coupled system (??). 
While / is maintained unchanged between both simulations, the distribution g does change. 
Thus, the corresponding H = f * g changes as it can be seen in part b; the blue, red, and 
green dots correspond to the speed change induced in an oscillator within the blue cluster 
by oscillators of each cluster. Since all clusters have the same number of oscillators, the net 
effect is zero. In part c the time evolution of oscillators' phases relative to the phase of a 
blue cluster oscillator are shown. Although the initial conditions are exactly the same, the 



wider delay distribution on Figure [13] produces negative slope on the red and green points 
of part b, which destabilizes the clusters and drives oscillators toward in-phase synchrony. 



Finally, we simulate the same scenario as in Figures 12 and 13 but now changing N and 
the standard deviation, i.e. the delay distribution width. Figure 14 shows the computation 
of the synchronization probability vs. standard deviation. The dashed line denotes the 
minimum value that destabilizes the equivalent system. As N grows, the distribution shape 
becomes closer to a step, which is the expected shape in the limit. It is quite surprising that 
as soon as the equilibrium is within the region of H with negative slope, the equilibrium 
becomes unstable as the theory predicts. 
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FIG. 12. Pulse-coupled oscillators with delay: Stable equilibrium 



V. CONCLUSION 



This paper analyzes the dynamics of identical weakly coupled oscillators while relaxing 
several classical assumptions on coupling, delay and topology. Our results provide global 
synchronization guarantee for a wide range of scenarios. There are many directions that 
can be taken to further this study. For example, for different topologies, to guarantee global 
in-phase synchronization, how does the requirement on coupling functions change? Another 
specific question is to complete the proof in Section HI C for the cases when m = 1, 3, 5. 
Finally, it will be of great interest if we can apply results and techniques in this paper to a 
wide range of applications such as transient stability analysis of power networks and clock 
synchronization of computer networks. 

Acknowledgments: The authors thank Steven H. Strogatz of Cornell for useful discussions. 



Appendix A: Proof of Proposition [3] 

As in Proposition [2] we will use our cut condition to show the instability of <fi* . Thus, 
we define a partition P = (S,V(G)\S) of V(G) by taking S to be a maximal subset of 
V{G) such that d{4>, S) < ^, see Figured for an illustration of P. Notice that any of these 
partitions will include all the oscillators of two consecutive blocks of every constellation. 

Instead of evaluating the total sum of the weights in the cut we will show that the sum 
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FIG. 13. Pulse-coupled oscillators with delay: Unstable equilibrium 



of edge weights of the links connecting the nodes of one constellation in S with the nodes 
of a possibly different constellation in V(G)\S is negative. In other words, we will focus on 
showing 

£ /U*-^)<0 (Al) 
where JC hh = {ij :ieC h nS,je C h n V(G)\S}. 



Given any subset of integers J, we define 

9 J m (S) = g m (S)-J2f(--j + S). 
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FIG. 14. Pulse-coupled oscillators with delay: Synchronization probability 



Then, we can rewrite (|A1|) as 



= (^°' 1} ) / (^y 2 ) + (^Ww,) 

27T 

=2^(5 W2 )-f(^ 2 + -)-2/ / (^ 2 ) 

lib 

- m lh - -) (A2) 

where Si t i 2 G [0, — ] is the phase shift between the two constellations. Then, if we can show 
that for all 5 G [0, — ] (A2) is less than zero then for any values of l\ and I2 we will have 
(Al ) satisfied. 

Since / is odd and even around |, /' is even and odd around | and g' m (S) can be rewritten 
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as 



+ E \nS+-j)-nS-8gn{j)- + -j)) 



i<lil<L|J 

+ —k) + f'(S — —k) 
m m 



-[k odd] 



where Q dd] is the indicator function of the event [A; odd], the sum is over all the integers 
j with 1 < \j\ < L|J and k = ^ 

The last term only appears when k is odd and in fact it is easy to show that it is always 
negative as the following calculation shows: 

-f{S+-k)-f{S--k) = 
m m 

= -f^k + 5)-f'(-k-5) 
m m 

= f(o)-f(e-<i>)<o 

where in step one we used the fact of /' being even, in step two we used k = and in step 
three we use /' being odd around |. The last step comes from substituting 6 = | — 5 + 
= ^ and apply Claim [3j since for m > 7 we have < 8 — <fi < 8 < ir. 

Then it remains the show that the terms of the form f'(5 + —j) — f'(5 — sgn(j)^ + —j) 
are negative for all j s.t. 1 < \j\ < [fj- This is indeed true when j is positive since for all 
Se [0,f] we get 

it 2n ^ 2n „ . k . 

< 5 1 j < 5 H j < 7r, for 1 < j < I -J 

m m m 2 

and thus we can apply again Claim [3j 

When j is negative there is one exception in which Claim [3] cannot be used since 

-7T < 5 + — j < 5 + — j + - < oys E [0, — ] 

m m m m 

only holds for — |_|J < j < —2. Thus the term corresponding to j = —1 cannot be directly 
eliminated. 

Then, by keeping only the terms of the sum with j — ±1, g' m is strictly upper bounded 
for all 6 e [0, f] by 

<4W </'(5) + - ^) - /'(^ - ^) 

+ f\5+-)-f'(5+-) (A3) 
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Now substituting ( |A3[ ) in (A2) we get 

E fw-M 

< /'(* - -) - 2f'(5 --) + f(5 +-)- 2f'(5 + *) 
m m m m 

<f{6-—)-2f'{6--)-f'(6 + -) 
m mm 

< f{8 -—)- 2f(5 - -) 
m m 

where in the last step we used the fact that for m > 6 and 5 G [0, ^], f{5 + > 0. 

Finally, since for 5 G [0, — ] f'(8 — — ) is strictly increasing and f'(S — — ) achieves its 
minimum for 5 G {0, then 

f(S --)- 2f'(5 --)< f(0) - 2f(-) < 
mm m 

where the last inequality follow from Claim |4j 

Therefore, for all m odd greater or equal to 7 we obtain 



E fty - ft) < /'(°) - 2 /'(-) ^ 

«jSAC i:l ; 2 



and since this result is independent on the indices Zi, I2, then 

E /'(^ - ft) 

ijeC(s,v(G)\s) 

= EE E m-ti)<° 

and thus 6* is unstable. 



[1] A. T. Winfree, "Biological rhythms and the behavior of populations of coupled oscillators," 

Journal of Theoretical Biology, vol. 16, no. 1, pp. 15 - 42, 1967. 
[2] C. S. Peskin, Mathematical aspects of heart physiology. New York, NY, USA: Courant Institute 

of Mathematical Sciences, New York University, 1975. 
[3] P. Achermann and H. Kunz, "Modeling circadian rhythm generation in the suprachiasmatic 

nucleus with locally coupled self-sustained oscillators: Phase shifts and phase response curves," 

Journal of Biological Rhythms, vol. 14, no. 6, pp. 460 - 468, 1999. 
[4] J. Garcia-Ojalvo, M. B. Elowitz, S. H. Strogatz, and C. S. Peskin, "Modeling a synthetic 

multicellular clock: Repressilators coupled by quorum sensing," Proceedings of the National 

Academy of Sciences of the United States of America, vol. 101, no. 30, pp. 10955-10960, 2004. 
[5] S. Yamaguchi, H. Isejima, T. Matsuo, R. Okura, K. Yagita, M. Kobayashi, and H. Okamura, 



26 



"Synchronization of cellular clocks in the suprachiasmatic nucleus," Science, vol. 32, pp. 1408- 
1412, Nov. 2003. 

[6] I. Z. Kiss, Y. Zhai, and J. L. Hudson, "Emerging coherence in a population of chemical 
oscillators," Science, vol. 296, pp. 1676-1678, Nov. 2002. 

[7] R. York and R. Compton, "Quasi-optical power combining using mutually synchronized os- 
cillator arrays," IEEE Transactions on Microwave Theory and Techniques, vol. 39, pp. 1000- 
1009, Jun. 1991. 

[8] Y.-W. Hong and A. Scaglione, "A scalable synchronization protocol for large scale sensor 
networks and its applications," IEEE Journal on Selected Areas in Communications, vol. 23, 
pp. 1085-1099, May. 2005. 

[9] G. Werner-Allen, G. Tewari, A. Patel, M. Welsh, and R. Nagpal, "Firefly-inspired sensor net- 
work synchronicity with realistic radio effects," in SenSys: Proceedings of the 3rd International 
Conference on Embedded Networked Sensor Systems, (New York, NY, USA), pp. 142-153, 
ACM, 2005. 

[10] S. A. Marvel and S. H. Strogatz, "Invariant submanifold for series arrays of josephson junc- 
tions," Chaos, vol. 19, p. 013132, Mar. 2009. 

[11] P. C. Bressloff and S. Coombes, "Travelling waves in chains of pulse-coupled integrate-and-fire 
oscillators with distributed delays," Phys. D, vol. 130, no. 3-4, pp. 232-254, 1999. 

[12] P. Ashwin and J. W. Swift, "The dynamics of n weakly coupled identical oscillators," J. 
Nonlinear Sci, vol. 2, no. 1, pp. 69-108, 1992. 

[13] G. B. Ermentrout, "Stable periodic solutions to discrete and continuum arrays of weakly 
coupled nonlinear oscillators," SIAM J. Appl. Math., vol. 52, no. 6, pp. 1665-1687, 1992. 

[14] Y. Kuramoto, "International symposium on mathematical problems in theoretical physics," 
in Lecture notes in Physics, vol. 39, p. 420, Springer, 1975. 

[15] O. V. Popovych, Y. L. Maistrenko, and P. A. Tass, "Phase chaos in coupled oscillators," Phys. 
Rev. E, vol. 71, p. 065201, Jun. 2005. 

[16] E. Brown, P. Holmes, and J. Moehlis, "Globally coupled oscillator networks," in Perspec- 
tives and Problems in Nonlinear Science: A Celebratory Volume in Honor of Larry Sirovich, 
pp. 183-215, Springer, 2003. 

[17] A. Jadbabaie, N. Motee, and M. Barahona, "On the stability of the kuramoto model of 
coupled nonlinear oscillators," in Proceedings of the American Control Conference., vol. 5, 
pp. 4296-4301, June 30, 2004. 

[18] D. Lucarelli and I. -J. Wang, "Decentralized synchronization protocols with nearest neighbor 
communication," in Proceedings of the 2nd International Conference on Embedded Networked 
Sensor Systems, 2004. 

[19] P. Monzon and F. Paganini, "Global considerations on the kuramoto model of sinusoidally 
coupled oscillators," in Proceedings of the 44th IEEE Conference on Decision and Control, 
and European Control Conference, (Sevilla, Spain), pp. 3923-3928, Dec. 2005. 

[20] R. E. Mirollo and S. H. Strogatz, "Synchronization of pulse-coupled biological oscillators.," 
SIAM J. Appl. Math, vol. 50, pp. 1645-1662, 1990. 

[21] L. Moreau, "Stability of continuous-time distributed consensus algorithms," in Decision and 
Control, 2004. CDC. 43rd IEEE Conference on, vol. 4, pp. 3998 - 4003 Vol.4, dec. 2004. 

[22] W. Ren and R. Beard, "Consensus of information under dynamically changing interaction 
topologies," in American Control Conference, 2004- Proceedings of the 2004, v °l- 6, pp. 4939 
-4944 vol.6, 30 2004-july 2 2004. 

[23] G.-B. Stan and R. Sepulchre, "Dissipativity characterization of a class of oscillators and net- 



27 



works of oscillators," in Decision and Control, 2003. Proceedings. 42nd IEEE Conference on, 

vol. 4, pp. 4169 - 4173 vol.4, dec. 2003. 
[24] L. Scardovi, A. Sarlette, and R. Sepulchre, "Synchronization and balancing on the n-torus," 

Systems & Control Letters, vol. 56, no. 5, pp. 335 - 341, 2007. 
[25] R. Sepulchre, D. A. Paley, and N. E. Leonard, "Stabilization of planar collective motion with 

limited communication," Automatic Control, IEEE Transactions on, vol. 53, no. 3, pp. 706- 

719, 2008. 

[26] A. Papachristodoulou and A. Jadbabaie, "Synchonization in oscillator networks with heteroge- 
neous delays, switching topologies and nonlinear dynamics," in Proceedings of the 45th IEEE 
Conference on Decision and Control, pp. 4307-4312, Dec. 2006. 

[27] E. M. Izhikevich, "Phase models with explicit time delays," Phys. Rev. E, vol. 58, pp. 905-908, 
Jul. 1998. 

[28] E. M. Izhikevich, "Weakly pulse-coupled oscillators, fm interactions, synchronization, and 
oscillatory associative memory," IEEE Transactions on Neural Networks, vol. 10, pp. 508- 
526, May. 1999. 

[29] E. Mallada and A. Tang, "Synchronization of phase-coupled oscillators with arbitrary topol- 
ogy," in Proceedings of American Control Conference, 2010. 

[30] E. Mallada and A. Tang, "Weakly pulse-coupled oscillators: Heterogeneous delays lead to 
homogeneous phase," in IEEE Conference on Decision and Control, 2010. 

[31] U. Ernst, K. Pawelzik, and T. Geisel, "Synchronization induced by temporal delays in pulse- 
coupled oscillators," Phys. Rev. Lett., vol. 74, pp. 1570-1573, Feb. 1995. 

[32] U. Ernst, K. Pawelzik, and T. Geisel, "Delay-induced multistable synchronization of biological 
oscillators," Phys. Rev. E, vol. 57, pp. 2150-2162, Feb. 1998. 

[33] Y. Kuramoto, Chemical Oscillations, Waves, and Turbulence. Berin Heidelberg New York 
Tokyo: Springer- Verlag, 1984. 

[34] R. Sepulchre, D. A. Paley, and N. E. Leonard, "Stabilization of planar collective motion: All- 
to-all communication," Automatic Control, IEEE Transactions on, vol. 52, no. 5, pp. 811-824, 
2007. 

[35] H. Daido, "Quasientrainment and slow relaxation in a population of oscillators with random 
and frustrated interactions," Phys. Rev. Lett, vol. 68, pp. 1073-1076, Feb 1992. 

[36] A. Rantzer, "A dual to lyapunov's stability theorem," Systems and Control Letters, vol. 42, 
pp. 161-168, 2001. 

[37] B. Bollobas, Modern Graph Theory. New York: Springer, 1998. 

[38] C. Godsil and G. Royle, Algebraic Graph Theory. New York: Springer, 2001. 

[39] H. K. Khalil, Nonlinear systems; 3rd ed. Prentice-Hall, 1996. 

[40] F. C. Hoppensteadt and E. M. Izhikevich, "Weakly connected neural networks," 1997. 
[41] A. Sarlette, Geometry and symmetries in coordination control. PhD thesis, University of Liege, 
Belgium, January 2009. 

[42] C. van Vreeswijk, L. Abbott, and G. B. Ermentrout, "When inhibition not excitation syn- 
chronizes neural firing," Journal of Computational Neuroscience, vol. 1, no. 4, pp. 313-321, 
1994. 

[43] W. Gerstner, "Rapid phase locking in systems of pulse-coupled oscillators with delays," Phys. 

Rev. Lett, vol. 76, pp. 1755-1758, Mar 1996. 
[44] S. Watanabe and S. H. Strogatz, "Constants of motion for superconducting josephson arrays," 

Phys. D, vol. 74, no. 3-4, pp. 197-253, 1994. 



28 



